Numerical comparison of different algorithms for construction of wavelet matrices 
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Abstract. Factorization of compact wavelet matrices into primitive ones has been known for 
more than 20 years. This method makes it possible to generate wavelet matrix coefficients 
and also to specify them by their first row. Recently, a new parametrization of compact 
wavelet matrices of the same order and degree has been introduced by the last author. This 
method also enables us to fulfill the above mentioned tasks of matrix constructions. In 
the present paper, we briefly describe the corresponding algorithms based on two different 
methods, and numerically compare their performance. 
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1. Introduction 



An m x (N + l)m matrix 
(1) A={A A 1 
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(Aj are square blocks) is called a wavelet matrix [6] if it satisfies the so called shifted 
orthogonality condition: 



5koI m , k — 0,1, 



N-k 

(2) E 

3=0 

where A* denotes the conjugate transpose of A, o" fco is the Kronecker delta, and I m is 
the m x m unit matrix. 

In the polyphase representation of matrix A, 

N 

(3) A(z) = Y,A k z k =:{^(z)}™ =1 , 

k=0 

the condition (T5]) is equivalent to 

(4) A(z)A(z) = I m , 

where A(z) = J2k=o A% z ^ k is the adjoint to A(z). 

In the sequel, the matrices of the form (pQ) and their polyphase representation ([3]) 
will be identified. 

Our notion of a wavelet matrix is weaker than usual as some linear condition is also 
required to be satisfied (see, e.g. [7]) which is irrelevant in our consideration. Instead, 

l 



2 



we require the condition 

(5) A(l) = / m . 

The integers m and iV are called, respectively, the rank and the order of a wavelet 
matrix (JTJ or ([3]) (it is assumed that ^ 0). It follows from that detA(z) 
has always the form cz d , d > 0, |c| = 1, and the integer d is called the degree of 
A. The class of wavelet matrices of rank m, order N and degree d will be denoted 
by W(m, N, d). In addition, Wo(m, N, d) denotes the class of those A G W(m, N, d) 
for which fl5} holds, and Wi(m, A 7 ", d) denotes the class of those ^4. G Wo(m, A 7 ", d) for 
which the last row of A^ differs from the zero vector of C m . 

It can be proved that the degree of any wavelet matrix is grater than or equal to 
its order, i.e. d > N (see, e.g., (2J Lemma 1]) and d = N holds except for some 
degenerated cases (see [7J p. 58]). We call the case d > N singular as the uniqueness 
of solutions, which we are going to construct numerically, fails to hold in this situation 
[6j. Namely, we consider wavelet matrices from the class W(m, N, N). It differs from 
Wi(m, N, N) by unitary multipliers on the left and on the right (see [2]). 

A wavelet matrix ~V(z) of order and degree 1 is called primitive. It can be shown 
(see, e.g. [TJ p. 59], [5J, [2]) that every V(z) G W (m, 1, 1) has the form 

V(z) = I m - v*v - v*vz, 

where v = (v\, V2, ■ ■ ■ , v m ) G C m is a vector of the unit norm, vv* = 1. 

The following wavelet matrix factorization theorem was first proved in the yearly 
90's in a related theory of multirate filter banks [S]. We formulate it for nonsingular 
matrices 

Theorem 1.1. For any A(z) G Wo(m, N, N), there exists a unique factorization 

N 

i=i 

where each ~Vj(z) G W (m, 1,1). 

This theorem provides a possibility to generate wavelet matrices of arbitrary order. 
The computational complexity of this method and its numerical tests are described 
in the next sections. Theorem 1.1. helps also to solve the following wavelet matrix 
completion problem [3], [BJ: Given the first row of a wavelet matrix, find its remaining 
rows, i.e. if the first row of ([T]) is given which satisfies the shifted orthogonality 
condition 

(N+l-k)m 

(6) a j^j+km = S k0 , fe = 0,l,...,JV, 

3=1 

then one should find the remaining entries of A which results in wavelet matrix. We 
emphasize that this problem has a unique solution (up to certain unitary matrix) if 
we search A in W(m, N, N) (see [7J Th. 4.17]). In the next sections, we describe and 
test numerically the existing algorithm of such construction. 
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A new parametrization of nonsingular compact wavelet matrices appeared in [2] 
in the form of Theorem 1.2 below which gives a one-to-one continuous map between 
C N (™~ 1 ) fmdWi(m,N,N) 

Let P N := {El 

=0 CkZ k : Co, C\, . . . , Cjv G C} be the set of polynomials and P N :— 
{ J2k=i c k z ~ k '■ c ii c 2, ■ ■ ■ , cn £ C} (note that P N fl P N = {0} according to our nota- 
tion). Up(z) = J2k=-N°k zk , then [p(z)] _ = Efc=-w c ^ fc and + = Efc=o c fc z& - 
Theorem 1.2. Let A > 1. For any polynomials 

(7) e^, i = l,2,...,m-l, 
there exists a unique 

(8) A(z) G Wi(m, N, N) 
such that 

(9) C, l (z)a lj {z) + C,2{z)a 2j {z) + . . . + ( m _ 1 (z)a m _ 1J (z) + a^ j (z) G P^, j = 1,2, . . .,m. 

Conversely, for each A(z) satisfying (jSj), there exists a unique (m— l)-tuple of Laurent 
polynomials (J7|) such that OH]) holds. 

Further refinement of Theorem 1.2 enables us to solve the wavelet matrix completion 
problem as well (21 §5]. The exact formulas of these constructions and numerical tests 
of corresponding algorithms are given in Sections 3 and 4. 

In conclusion we analyze numerical performances of described algorithms and, based 
on these data, compare two different methods. 

2. The existing algorithms of wavelet matrix construction 

The following wavelet matrix generation procedure is based on Theorem 1.1 

Algorithm 2.1. Step 1. Take arbitrary nonzero vectors Vj- G C m , j = 1, 2, . . . , N, 
(they can be selected randomly) and let 

(10) i) (v /V ;) 'vjv,. 

Then ~Vj(z) — I m — Pj + Pjz, j = 1, 2, . . . , N, are primitive wavelet matrices. 
Step 2. Let A Q (z) = I m and for j = 1, 2, . . . , N do: 

(11) &j{z)=A j - X {z){l n -P j {I m -zj). 

(Matrix multiplication in ( II ip requires approximately m 2 (j — 1) operations (ops) count- 
ing only multiplications. Thus the cycle in Step 2 needs ~ Ej=i m2 {j ' ~ 1) = 0(m 2 N 2 ) 
ops.) 

Then A(z) = A^(z) will be the wavelet matrix of rank m and degree N . It can 
be seen that ord(^4) = N if and only if Vj-v* +1 ^ for j = 1,2, ... ,N — 1 (i.e. 
the consecutive Vj-s in (flO]) are not orthogonal), and in this case the last row of 
An = Wj=x Pj differs from G C m if and only if the last coordinate of vi differs from 
0. Thus, for randomly selected v^-s in Step 1, the wavelet matrix A(z) belongs to 
Wi(m, A, A) with probability 1. 

The following procedures describe a numerical solution to the wavelet matrix com- 
pletion problem [B], [3]- 



Algorithm 2.2. Given 

(12) a= (a[,al, ■■■ ,a\ N+i)m ) =: (ao,ai...,a w ), a N ^ 

satisfying conditions (jSJ) and ^^L a, = e\ — (1, 0, . . . , 0) G C m . 

Step 1. Let Fat = (a A ra^ r )~ 1 a^ r a A r and let (a^, . . . , a^) := (a , a x . . . , a^). 
For j = N, N - 1, . . . ,2 do: 

a a-D = gC,-) + (a w _ a a) )jR?5 i = o, i, . . . , iv — l, 

and 

r 3 ^\ — ^a 7V _ 1 ^a iV _ 1 ) ) \<*n-i ) d jv-i ■ 
(This step needs approximately 0{mN 2 ) ops.) 
Step 2. Compute the product 

N 

aw = n( j ™ - p i + 

i=i 

using Step 2 of Algorithm 2.1. 

Then ^4 = A(z) is the unique wavelet matrix from yVo(m, N, N) with the first row 
(fT2|) (see also [7, Th. 4.17]). 

All in all, the number of operations in Algorithms 2.1 and 2.2 can be estimated as 
0{m 2 N 2 ). 

3. New algorithms of wavelet matrix construction 

In this section we describe algorithms based on recently developed method of wavelet 
matrix parametrization [I]. First we generate A € Wi(m, N, N) (see [2] for justifica- 
tion of the given procedures). 

Algorithm 3.1. Step 1. Take arbitrary m — 1 Laurent polynomials from 

N 

(13) (i{z) = ^~f ik Z- k , l = 1,2,... ,771-1, 

fe=l 

(the coefficients 7^ can be selected randomly). 

Step 2. Perform upper triangular, diagonal, lower triangular factorization 

(14) A = UDU* 
of 

m—1 

(15) A = ^e J e~ + / JV+1 , 

i=l 

where Gj is the upper triangular (N + 1) x (N + 1) Hankel matrix with the first row 
(0,7ii>7i2, • • • ,7iJv)- 

Since A has a displacement structure of rank m (see [H Appendix]) the factorization 
(fll]) can be performed in 0{mN 2 ) ops (as it is described in [5j Appendix F.l] ) without 
constructing ( ITS]) explicitly 
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Step 3. Solve the system of (N + 1) x (N + 1) linear algebraic equations 
(16) AX = Bj 

to times taking different right hand sides, where Bj = (0, jn, 7^, • • • , JiN) T , j = 
1, 2, . . . , to - 1, and B m = (1, 0, . . . , 0) T 

Since we have the factorization (JHJ), the solution of the system ([16]) requires 0(N 2 ) 
ops and Step 3 totally needs 0(mN 2 ) ops. 

Let (ofjo, CKji, . . . , OLjN) be the solution of (fTBT) and let 
v 

U i( 2 ) = ft i fc2;_fc alld h rnj(z) = Z N VLj(z), j = 1, 2, . . . , TO. 
k=0 

Step 4- Compute the coefficients of the following polynomials from 

biji z ) = [G( z ) u j( z )] + - Sij, 1 <i <m, 1 < j < m. 

As the multiplication of polynomials of order N takes 0(N\ogN) ops by FFT, Step 
4 totally needs 0(m 2 N log N) ops. 

Step 5. Constructing the matrix polynomial B(z) = {by (z)}™j =1 , 

A(z)=B(z){B(l)y 1 

will be a wavelet matrix from Wi(to, N, N). 

Since m x to matrix inversion needs 0(m 3 ) ops, Step 5 totally needs 0(Nm 3 ) ops. 

Now we describe a new algorithm of wavelet matrix completion based on Theorem 
1.2. Its justification can be found in [2]. 

Algorithm 3.2. Data is the same as in Algorithm 2.2. 

Step 1. Select a coordinate of a N = (al nN+1 , Omiv+2' ' ' ' i a ln(N+i)) w ith maximum 
absolute value. Since a^r ^ 0, this coordinate differs from and let it be a^jv+j- This 
preliminary step will improve the accuracy of the final result. 

Step 2. Let (an(z), ai 2 (z), . . . , ai m (z)l be the polyphase representation of (TP2"]) . i.e. 
the first row of ()3]). 

Compute the first iV + 1 coefficients, say 70, 71, ... , 7tv, of the reciprocal (in a neigh- 
borhood of 0) of Yl!k=o^mk+j zN ~ k = zN ^ij{ z )-> where j was determined in Step 1. 
(This step requires 0(N 2 ) ops, though some papers pQ report that it can be done in 
0(N log N) ops using parallel computations.) 

Let C(z) = Y.k=nlkZ k . _ 

Step 3. Compute Q(z) = [a.u(z)£(z)] for j 7^ i = 1,2,..., to. (This step needs 
0(mN log N) ops.) 

Step 4- Use Algorithm 3.1 with the data (i(z), ( 2 (z), . . . , Q_i(z), Q + i(z), . . . , Cm( z ) 
to construct the corresponding wavelet matrix. Denote this matrix by Aj(z) (in 
polyphase representation). Then, if we transpose A}(z) and move its last row in 
the place of zth row and its last column in the place of ith column, we get the desired 
A(z) G Wo(to, N, N). 

All in all, the number of operations in Algorithms 3.1 and 3.2 can be estimated as 
0(TO 2 AnogAO +0(to 3 A). 
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4. Numerical Simulations 

To compare the performance of the described algorithms, their computer code was 
written in Mathematica-8. A PC with 2.40GHz Intel Quad Core CPU and 2GB RAM 
was used for numerical simulations. 

The accuracy level of the wavelet matrix construction algorithms (2.1 and 3.1) is 
measured by that of relation (4), in which we substitute the computation results, 
whereas the accuracy level of wavelet matrix completion algorithms (2.2 and 3.2) is 
naturally measured by the difference between the initial data and the first row of the 
computed matrix. 

The accuracy levels determined in this way are essentially the same for both methods 
and are quite close to precisions in which Mathematica-8 carries out calculations. (As 
it is known Mathematica-8 provides an opportunity to make this precision arbitrarily 
large.) Therefore, all experiments for comparison of the speeds of different algorithms 
were run in the standard double precision. The results of these simulations are pre- 
sented in the Tables below. As it was expected (since 0(m 2 N 2 ) should be larger than 
0(m 2 N log N) + 0(m 3 N) for m < N), the new algorithms are faster than the old 
ones, and the difference between their performance times is becoming more and more 
evident as m and N grow, keeping N sufficiently larger than m as it should be for 
wavelet matrices applicable in practice. 

Table I 

Results of Computer Simulations of Wavelet Matrix Construction Algorithms 



Rank m 


10 


10 


10 


10 


20 


20 


20 


30 


30 


30 


50 


50 


50 


Order N 


50 


150 


300 


400 


100 


250 


300 


100 


150 


200 


100 


200 


300 


Time (New Alg.) 


0.34 


2.80 


10.60 


18.86 


2.58 


14.93 


20.12 


4.52 


9.14 


15.75 


8.65 


27.93 


57.54 


Time (Old Alg.) 


0.79 


5.45 


16.09 


19.44 


10.84 


57.44 


62.52 


23.85 


49.98 


85.93 


67.69 


234.48 


401.16 



Table II 

Results of Computer Simulations of Wavelet Matrix Completion Algorithms 



Rank m 


10 


10 


10 


10 


20 


20 


20 


30 


30 


30 


50 


50 


50 


Order N 


50 


150 


300 


400 


100 


250 


300 


100 


150 


200 


100 


200 


300 


Time (New Alg.) 


0.39 


2.83 


10.96 


19.76 


2.85 


15.76 


22.77 


4.58 


9.46 


17.13 


9.23 


28.01 


58.48 


Time (Old Alg.) 


0.85 


6.08 


18.24 


22.01 


11.96 


63.13 


70.54 


24.96 


52.91 


88.09 


74.42 


333.71 


473.61 



5. Conclusion 



In this paper, we describe in detail two new algorithms of wavelet matrix construc- 
tion and completion, introduced in [2J. The results of numerical simulations are pre- 
sented, which prove the advantage of the new algorithms over the existing algorithms 
in performance speed. 

Another advantage of the new method, which should be mentioned here and might 
be used in the future, is that the algorithms based on this method can be divided into 
m parallel tasks which will make them even more faster. The old method is heavily 
recurrent and misses any opportunity to be parallelized. 
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